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Abstract 

In recent work we have shown how an accurate reduced model can be utilized 
to perform mesh refinement in random space. That work relied on the explicit 
knowledge of an accurate reduced model which is used to monitor the transfer 
of activity from the large to the small scales of the solution. Since this is 
not always available, we present in the current work a framework which shares 
the merits and basic idea of the previous approach but does not require an 
explicit knowledge of a reduced model. Moreover, the current framework can be 
applied for refinement in both random and physical space. In this manuscript we 
focus on the application to random space mesh refinement. We study examples 
of increasing difficulty (from ordinary to partial differential equations) which 
demonstrate the efficiency and versatility of our approach. We also provide 
some results from the application of the new framework to physical space mesh 
refinement. 

Keywords: Adaptive mesh refinement, Multi-element, gPC, uncertainty 
quantification, discontinuities. 


1. Introduction 

The prediction of the evolution of most physical systems is inevitably affected 
by uncertainties in the initial, boundary and/ or parameters entering the formu¬ 
lation of the governing equations. The two prevalent methods for quantifying the 
effect of this uncertainty on the evolution of the statistics are Monte Carlo (MC) 
and stochastic spectral method (see, e.g., [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] and 
references therein). Both of these methods lead to efficient algorithms as long as 
the number of sources of uncertainty (uncertainty dimensionality) is relatively 
small. However, for most realistic systems, the dimensionality of uncertainty is 
rather large. This can lead to very serious obstacles in the application of both 
MC and stochastic spectral method. Development of sparse grid techniques 
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[13, 2, 14, 15] as well as sparse polynomial bases [8] for the Galerkin formula¬ 
tion have greatly alleviated this difficulty, however these techniques heavy rely 
on the regularity of the solution in the random space. At the same time, not 
all values in the range of the sources of uncertainty contribute equally to the 
determination of the solution statistics. This creates the hope that if one can 
allocate more computational resources around the most important values of the 
random variables, then one can still produce reasonably accurate predictions 
of the solution statistics. Multi-element generalized polynomial chaos method 
(ME-gPC) [16, 17] multi-element probabilistic collocation method (ME-PCM) 
[18] and multi-resolution wavelet expansion [19, 20] provide the framework to 
obtain the statistics of the system when applying different computational re¬ 
sources to different area (element) in random space. In order to affect the 
concentration of resources in the more sensitive areas of the random space (for 
example, the discontinuities in the random space), one needs tools for the re¬ 
liable location of these sensitive areas. In particular, if we consider a mesh 
discretization of the random space, then one needs to be able to decide when, 
where and in which direction to refine this mesh. [16] provides an approach to 
refine the random elements and allocate the computational resources adaptively 
as the system evolves. Adaptive sparse grid methods [21, 22, 23] also provided a 
dynamical way to obtain unconstructed collocation mesh to resolve the system. 

In [24] we presented a novel way of performing mesh refinement in random 
space. The approach was based on the use of a reduced model of the full order 
model, when such a reduced model is available in explicit form. The main idea 
was that successful mesh refinement depends on the accurate monitoring of the 
transfer of activity (e.g. mass, energy) from the large to the small scales of the 
solution. At the same time, this is exactly the role of a good reduced model, to 
reproduce faithfully this transfer of activity. Thus, if one has access to a good 
reduced model, then one can use it to perform the mesh refinement task. In [24] 
we provided theoretical results which established the soundness of our approach 
as well as numerical experiments which corroborated the theory. However, the 
whole framework relied on the explicit knowledge of the reduced model which 
is not always easy to obtain. 

In the current work, we present a simple reformulation of the idea in [24] 
which retains the merits of the previous approach while at the same time elimi¬ 
nating the need for the explicit knowledge of the reduced model. In particular, 
we offer a way to compute the rate of transfer of activity from the large to the 
small scales by using appropriate quantities which can be computed without the 
explicit knowledge of the reduced model (see Section 3). The main idea is that 
at every step and at every element where refinement may be needed, one can 
obtain an expansion of the solution (in the element) on some chosen basis. The 
coefficients of this expansion can be used to compute the rate of transfer of 
activity from the large to the small scales (again within the element) and thus, 
decide whether there exists a need for refinement. We present results for the 
uncertainty quantification of three different cases of increasing complexity (lin¬ 
ear ODE with uncertain coefficient, Kraichnan-Orszag system and Kuramoto 
Sivashinsky equation with an uncertain bifurcation parameter). In all three 
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cases, the mesh refinement algorithm is capable of detecting the sensitive areas 
of the random space with high accuracy. 

The adaptive mesh refinement framework described above can also be ap¬ 
plied to situations where refinement is required in physical space. The only 
difference with the random space case is that in physical space when assigning 
values for the modes in the newly generated elements, the coupling between 
adjacent elements is required. In other words, the elements are not decoupled 
and the continuity (or more) of the solution across the boundary of adjacent el¬ 
ements should be guaranteed. In the current work we present numerical results 
for the mesh refinement around the location of shocks for the ID Burgers equa¬ 
tion while more elaborate applications will appear in a forthcoming publication. 

2. Representation of uncertainty 

Let (£l,A,V) be a complete probability space, where D is the event space 
and V is the probability measure defined on A G 2 n , the a— algebra of subsets 
of D (these subsets are called events). Let D be a subset of R c (c G {1,2,3}) 
with boundary dD. Let C be an operator on D , and consider the following 
stochastic differential equation: 

«t(x, t;w) = £(x, t,u>\u), x G D. (1) 

The operator C usually involves spatial derivatives and can be either linear or 
nonlinear. In addition, C may depend on w G Cl. Appropriate initial conditions 
and boundary conditions sometimes involving random parameters relying on u 
are assumed so that the problem (1) is well-posed V—a.e. We also assume that 
for V—a.e. solution u := u(x, t\ui) is a function taking values in R. 

The random dependence operator C must satisfy a few important properties. 
The most important one is the ’’finite dimensional noise assumption” [9, 2], that 
is the random input can be represented with a finite-dimensional probability 
space. More precisely, the random input can be represented by a finite set of 
random variables, e.g., £(w) = (£i(w),--- ,£d(w)) be a d-dimensional random 
vector for ui G Cl. Then by the Doob-Dynkin lemma [25], the solution u(x, f;w) 
can be written as it(x, t; £(w)). 

Any second-order stochastic process can be represented as a random variable 
at each spatial and temporal location. Applying the Karhunen-Loeve expansion 
[6] along with the ” finite dimensional noise assumption”, a second-order stochas¬ 
tic process can be characterized by a finite set of mutually independent random 
variables. Thus without loss of generality, in this work we assume {£j(u;)}({_ 1 
are independent random variables. Let T := Ilf =1 T,, where T^ is the image of 
£i(f2), for i = 1,..., d. Let p(£) be the probability density function(p.d.f) of £. 
Then the problem (1) can be restated as following: find u : D x T —> R such 
that p-a.e. for £ G T, the equation 


«t(x,f;£) = £(x,t,£;u), x G D, 
with appropriate initial and boundary conditions. 


( 2 ) 
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2.1. Multi-element decomposition in random space 

To deal with irregularity (e.g. discontinuities) in the random space depen¬ 
dence we discretize the parametric space T into a collection of non-overlapping 
hypercubes. Since any d-dimensional independent random vector can be trans¬ 
formed to a d-dimensional independent uniform random vector, without loss of 
generality, let £ be a d-dimensional random vector defined on the random space 
T, where £j,i = 1 ■ • • ,d are independent identically distributed (i.i.d) uniform 
random variables defined on [—1,1] with p.d.f fi = Then T = [—1, l] d C 
with joint p.d.f. p = Following the notation in [16], we briey introduce the 
decomposition of the uniform random space (see [16] for more details). Let T 
be decomposed into TV non-overlapping hypercubes as following: 

Bk = [af, b\) x [a^, 62) x • • • x [a^, 6^], 

N 

r = U B k , (3) 

/c=1 

Bi n Bj = 0 if i / j, 


where i : j,k = 1,2,-•• ,7V. Let Xk,k = 1,2, ••• ,7V be the indicator random 
variables on each of the elements defined by 



if£eB fc , 

otherwise. 


For each random element B k , the local random vector is defined by 


C fc = (Ci\C 2 V--,C*) = £U, 


subject to the conditional p.d.f 


fc k 


1 

2 d Prob( Xfc = 1) ’ 


As = 1,2,---,7V, 


where Prob(xfc = 1) = 6 * 2 gi ■ After that, we transfer each to a new 

random vector £ fc defined on [—1, l] d by a map gk , 


9k(C k ) ■■ $ 


b k - a k k b\ + a k 
2 ^ + 2 


i — 1, 2, • • • , d. 


Thus, £ k is a random vector defined on [—1, l] d with constant p.d.f f k = 
With this decomposition, we can solve a system of differential equations with 
random input £ by combining the local approximations via subject to a 
conditional p.d.f. In practice if the solution u(£) is locally approximated by 
Uk(C k ),k = 1,2, ,7V, then the mth moment of u(£) on the entire random 
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space can be obtained by 


/*m(a(0) 


/ r »”(Op<i« 

N r 

^Prob( Xfc = l)/ uT(C k )fc*dC k 


fc=l 

AT 


E Prob (» = l)E(w m (x, t- -)|£ e B k ) 

k =1 

N r 

^prob( Xfe = i) / uT(9^(e))f k (e)de 
./r 


fc=i 


(4) 


U.i?. Localized stochastic discretization 

To solve the equation (1), any stochastic method (such as Monte Carlo sim- 
ulations(MCS), general polynomial chaos (gPC) Galerkin methods, stochastic 
collocation methods and etc.) can be applied in random space and result in a set 
of deterministic equations in the physical space which can be solved by standard 
numerical techniques. In the multi-element framework, MCS is straightforward, 
while both gPC Galerkin and stochastic collocation methods can be adapted to 
fit in. 

Next, we will briefly introduce the two most popular intrusive methods to 
deal with random inputs in one element. 

2.2.1. ME-gPC 

Without loss of generality, consider an orthonormal generalized polynomial 
chaos basis {<3>i}j??_ 0 spanning the space of second-order random processes on 
the probability space element B k (i = (*i, ■ ■ ■ ,id) G Nq is a multi-index with 

|i| = *i +-1- id)- The basis functions (£ fc (cn)) are polynomials of degree |i| 

with orthonormal relation 

(5) 

where dy is the Kronecker delta and the inner product between two functions 
g(£) and h(£) is defined by 

(g,h) Bk = [ g(C k MC k )fc-(C k )dC k • (6) 

J B k 

The solution u is approximated by the truncated gPC expansion 

U (x,t;C fe )= (7) 

|i|—0 

Substituting equation (7) into the governing system (1), we obtain the following 
system 

p ftf.k p 

*?)■ (8) 

|i|=o |i|=o 
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By using the Galerkin projection of (8) onto each element with the orthonormal 
polynomial basis { < hf}j , i | =0 , we derive 

f)v k p 

= |j| — 0,1, - - - ,p. (9) 

|i|—0 

This is a set of n = ( d ^ p ) coupled deterministic equations of the random modes 
ft; (x, t), |i| = 0,1, • • • ,p. Techniques for deterministic equations can be imple¬ 
mented to solve this system of equations. Then the mean of the solution on 
element B k can be evaluated by 

E(u(x,i; -)l£ G B k ) = f u(x,t;£))f(£)d$, 

JB k 

P (10) 

|i|=0 



The m-th moment of the approximated solution u(x, t: £) over element B k can 
be evaluated similarly. The global m-th moment of u(x, t;£) can be obtained 
by (4). 

2.2.2. ME-PCM 

To avoid solving a large coupled system of deterministic equations one can 
use the multi-element polynomial collocation method (ME-PCM) [18]. To im¬ 
plement ME-PCM, we need to specify the set of collocation points. Once the 
mesh discretization of T is prescribed, a set of collocation points {z k }j=i tr is 
identified in each element Bk, where r is the number of points used. Usually, 
these points are chosen to be the points of a cubature rule on Bk with integra¬ 
tion weights {w k }r_- l, k = 1 ,... ,N. Here we consider full tensor products of 
Gauss quadrature points for d = 1,2, and sparse grid points for d > 3. 

At each collocation point z k , we find the solution of the deterministic prob¬ 
lem 

u t (x.,t;z k ) = C(x,t,z k -u), x e D, (11) 

with appropriate initial and boundary conditions. Then for each element Bk, 
k = 1,..., N, we can construct an approximate solution I^ k u(x, t; £) using the 
set of solutions on the collocation points. For instance the operator can be 
chosen to be the tensor product Lagrangian interpolant, i.e., 

r 

I^ fe u(x,f;£) = ^2u(x,t-z k )£ k (C k ), (12) 

i =i 

where £ k (C k ) is the Lagrange polynomial corresponding to the collocation points 
{z k Y j=1 , refer [18] for details. Another frequently used operator is the orthogo¬ 
nal polynomial interpolant (gPC interpolant). To do this interpolation, assume 
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in element B k , u(x,t; £) can be approximated by u(x, t; £) = ^ m^(x, t)$ k (£ k ), 

|i|—0 

then by the cubature rule and the orthonomality of $*(C fc ), we have 

r 

«f(*,x) = ^ U (t,x; z*)$f(z*)u£, |i| = 0, ■ • ■ ,p. (13) 

i=i 

The gPC interpolation is constructed by 

4»«(x.*;€) = E (i4) 

|i|=0 


As a result, the global approximation is defined as: 

N 

u(x,i;£) = ^r. (15) 

fc=l 


where i = 1 refers to Lagrangian interpolation and i = 2 refers to gPC interpo¬ 
lation. 

Considering the statistical properties of the solution zt(x, t; z), once the m-th 
moment of u is calculated on each stochastic element the global m-th moment 
can be assembled by (4). The conditional m-th moment on element B k can be 
evaluated via the cubature rule over B k by 


E(u m (x,t;-)|^)eB t )= [ uZ(x,t;C k )fc*dC k 

■>B k 


E- m ( 


c i t; z k )i 


1=1 


(16) 


Again, we can obtain the global m-th moment by (4). For the properties of this 
discretization and error analysis please refer to [18]. 


3. Adaptive mesh refinement 

The structure of the division of the random space as described in the previous 
section can be allowed to change with time. In particular, different regions 
of random space can require different resolutions. This situation requires an 
adaptive mesh refinement procedure. In [24] we introduced an adaptive mesh 
refinement algorithm based on the reduced model of the full system. Once 
a reduced model is obtained, a quantity guiding the refinement is calculated. 
Unfortunately, a good reduced model of a general system is seldom available in 
an explicit form. However, the idea to use the memory terms in reduced model 
inspires the quantity used in this paper for the refinement. 

To convey the main idea, we can focus on one random space element. Sup¬ 
pose that the truncated Galerkin projection of equation (1) with orthonormal 
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basis is available, and assume p is the highest degree of the polynomial basis in 
the truncated system, then we obtain the full model 

dv- p 

= (£(x,i,£; ^2 Ui$i),$j), |j| = 0,1,--- ,p. (17) 

I i|—0 


For the reduced system, choose po < p (for our numerical examples we have 
chosen the ratio po/p to be around 1/2). 

Then the truncated system up to po can be written as: 

r)~ Po 

= (£(x,f,|; ^ ui$i),$j>, |j| = 0, l,--- ,p 0 . (18) 

|i|—0 


Let 

Po Po 

EpoW == E = II E WOllr, 

|j|=o I j I—o 

and 

Po Po 

E po (ti) := Y, Sj = II E 

IJI=o I j I—o 

E„ (u) can be interpreted as the L 2 norm of the projection of u onto the space 
spanned by {$j(£)}|j|=o- Then the rate of change of E po (ft) — E po (ft) denoted 
by Q is obtained as: 


n _^(E Po (m) — E po (u)) _ 0 ^ diij „ 0 ^ diij _ 

V- - dt u i 1 2^ dt u i 

|j|=o Li l=o 


dt 


(19) 


=2 E (A x ,t,£; E - 2 E (£( x >£>£; E wi$i),$j)fij. 

|j|=0 I i| =0 |j|=0 I i|—0 


Q measures the energy transfer between the full system (truncated up to order 
p) and the reduced system (truncated up to order p 0 ). Let Q = J D \Q\dx. When 
the quantity Q multiplied by the volume of a specific element exceeds a user- 
prescribed tolerance then the element is split in two parts i.e., mesh refinement 
takes place (see also [24]). 

If the random space has d dimensions with d > 2 then we need to decide 
not only when to refine but also in which directions. This can be decided in 
different ways and here we offer two possibilities. 

The first is to consider the contribution of the highest degree (po) basis in 
that dimension in the reduced model. In order to implement we define 

r P 

Sj = / |2(£(x, t, £; 'y ' Uid^i) , $p 0 ei )u po ei 
Jd |i|=0 

„ < 2 °) 

— 2(/2(x, t, {•; y ' {ti"!?;), $p oei )Cip oei | dx, 

|i|=0 
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where i refers to the z-th dimension, and is the index vector with z-th dimen¬ 
sion 1 others 0. Compared with (19), sj is the contribution in degree po basis 
in the z-th dimension. 

The second criterion for directional refinement is related to the total variance 
in the z-th direction. We define 


— I | ^ ^ ' (/l(x, t , ^ )&nei 

Jd n=l |i|—0 

Po P 

2^^{H(x, f, ^ ( Ui^j), )li n ei | dx. 

n=l |i|=0 


( 21 ) 


The quantity sf is the contribution to Q from all the degrees 1 to po in the z-th 
dimension. 

We now have all the ingredients needed to present our mesh refinement 
algorithm: 


Mesh refinement algorithm 

Step 1 Choose values TOL\ > 0 and TOL 2 > 0 for the tolerances. 
Step 2 Mesh refinement: 

For time step t <—!,■■ ■ , Nt 


Loop over all elements: 

On the fc-th element B update the modes for Bk 
If QPr(13fe) > TOLi, and if d = 1, split element B k and update 
system, otherwise loop over all dimensions: 

If si ( or sf, respectively) > TOL 2 -niax :)= i i ... ^ s]( or s|, respectively), 
split the element B & in two equal parts along the zth di¬ 
mension and update the system 
End if 
End if 

Update the information of the new elements 
End loop 

We have to make several remarks about the implementation details of the algo¬ 
rithm. 

First, in (19) we obtain an expression for the rate of change(<5) which does 
not involve temporal derivatives. There are two different ways to obtain Q 
which are based on the method used to solve (1). If the gPC Galerkin method 
is employed to solve (1), then at every time step the coefficients in the gPC 
expansion are available. In this case, Q is computed directly via (19). If the 
gPC collocation method is utilized to solve (1), then extra steps to evaluate the 
coefficients of the gPC expansion (see (13)) need to taken. However, since when 
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applying collocation method the resulting system is decoupled, the collocation 
approach is usually more efficient than the Galerkin approach. 

Secondly, whether one uses a Galerkin or a collocation method, there are 
two ways to compute Q then Q. If we are employing a Galerkin method, then 
the first way is to evolve, for each element Bk , both the full and reduced systems 
of equations i.e., solve equations (17) and (18). On the other hand, If we are 
employing a collocation method, then for each element Bk we use the tensor 
mesh with r points in each dimension and evolve the system on this mesh. 
Afterwards, we use a cubature rule to get the coefficients u\, |i| = 0,... ,p via 
(13). Here r depends on the highest order of polynomial representation of each 
random variable in the RHS of (9). To get the coefficients of the reduced system 
in element Bk, we use another tensor mesh with r o points in each dimension 
and evolve the system, then evaluate the coefficient u k for polynomial bases 
up to order po- For both approaches, we obtain Q by (19). If the explicit 
expressions for both full system and reduced systems (depending only on the 
coefficients of the expansion in the orthonormal basis) are available, then only 
p + 1 collocation points in each dimension for full system and po + 1 collocation 
points for the reduced system are required. If we are employing a Galerkin 
method, the second way to obtain the quantity Q is to evolve only the full 
system (9), and let u k = u\ , for |i| = 0,... ,po- If we are using a collocation 
method then we only evolve the system on a tensor mesh with r points in each 
dimension. Evaluate the coefficients for u k , |i| = 0,... ,p, then follow the steps 
as in the Galerkin case to get Q. 

The third remark concerns the generation of values in the newly created 
elements after splitting of the original element. For the Galerkin approach, 
when the element B k is split into two elements B kl and B k2 , two new random 
variables £ kl and 2 are generated according to the procedure in [16]. For 
the collocation approach, after the refinement a new mesh is generated. Let 
{zj 1 }][—i and {z^} r j=1 be the collocation points in the newly generated elements 
B kl and B k2 , respectively. We can use interpolation (see (12) or (14)) to obtain 
the system values u(x, r; z kl ) and u(x, r; z^' 2 ), j = 1,..., r (here r is the time 
when the refinement occurs). For the numerical examples below, we use gPC 
interpolation to update the system values on the new collocation points, since 
the gPC expansion is already available after evaluating Q. 

Our final remark is very important because the performance of the mesh 
refinement procedure relies heavily on it. For many cases, starting with only 
one element in random space is not enough, even if one uses high order gPC. 
If one insists on starting with one element, then the new points generated after 
the mesh refinement are not well placed to guarantee the accurate evolution of 
the solution. This is to be expected since an under-resolved representation of 
the random space initially will lead to erroneous identification of the regions 
that need to be refined. To avoid this one must have an adequate number of 
elements initially. If there is no other way to estimate how many elements are 
needed initially, one can decide by running several realizations with a different 
initial number of elements and monitoring the convergence of the results as the 
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initial number of elements is increased. 


4. Numerical Examples 

We have tested our adaptive mesh refinement algorithm both with a stochas¬ 
tic Galerkin method (AMR-GA) and with a stochastic collocation method (AMR- 
CO) on variety of problems of increasing complexity. In particular, we have 
applied it to a single linear ODE with a one-dimensional random parameter, 
the nonlinear Kraichnan-Orszag three-mode system with random initial con¬ 
ditions, as well as the ID Kuramoto-Sivashinsky (K-S) equation with random 
bifurcation parameter. 


4-1- One-dimensional linear ODE 

We begin by considering the simple ODE [16] 

^- = - k ( u )) u , u(0;w) = u 0 , 

at 

where k(uj) ~ U(— 1,1). The exact solution of this equation is 

u(t, ui) = uoe~ K '^ t . 


The statistical mean of the solution is 


Ha fp* _ p-t) f > n 

At(u(i;w)) = { 241 ’ 


u 0 , 


t = 0, 


and the variance is 


a 2 {u(t]u)) = { 4t V ) 3FD 


^^ 2t -Le- 2t - 2 ), t > 0, 


t = 0. 


( 22 ) 


(23) 


We will focus on the implementation of the adaptive mesh refinement with 
gPC collocation approach without an explicit expression of Galerkin projection. 
Also, we will evolve only the full system and use that to compute all the neces¬ 
sary quantities to decide when and where to refine. Since the random input is 
uniformly distributed, we choose Gaussian quadrature points as the collocation 
points. Let p be the order of the full system and po = \^-~\ the order of the 
reduced system. From Gaussian quadrature we know that if the gPC expansion 
of the solution is up to order p , then the highest polynomial order in random 
variable z of RHS in (26) for which the quadrature is exact equals 2p + 1. This 
means that in order to evaluate exactly the integral of these polynomials we 
need p+ 1 collocation points in each dimension of each element. Let be 

the collocation points in element B k , and {w k }jSf are the Gaussian quadrature 
weights respectively. Also, let u(t, z k ), i = 1,... ,p + 1 be the solutions to (22) 
at points f at time r. Let 

= (24) 

i =0 
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where {4>i} are the orthonormal Legendre polynomial basis functions. The 
coefficients uf(r), i = 1,... ,p. can be obtained through (13). To estimate Q, 
first we substitute (24) to left hand side(LHS) of (22) and we obtain: 


duf 


=- zu (t, z ), i = 0,---,P- (25) 

2=0 

By the orthonomality of i = 0,... ,p, we have the following expression: 

/ zu(t, z)$i(z)dz, i = 0,...,p. (26) 

J B k 


d*t 

dt 


By the quadrature rule, at time r, (26) becomes 


du? 

dt 


p +i 


= - E Z j U ( T ’ ) W j 1 * = 0, • • • ,P- 

3 =1 


(27) 


Immediately, we have: 


d_ 

dt 


E ("?) 2 

2—0 


2 E 



Po P+1 

- 2 EE *i z j u ( T ’ 

i=0 j=1 


(28) 


As mentioned in the previous paragraph we do not evolve the reduced truncated 
system up to order p 0 . Instead, we let v%(t) = Uj(r), for i = 0,... ,pq. To eval- 

Po 

uate the RHS of (18), we substitute the approximation u(r, z^) = 

2—0 

j = 1,... ,p + 1 and apply the quadrature rule for the integral. We find 


d 

dt 


Po 

E 

2—0 


R fc ) 2 



2=0 


po p+1 po 

- 2 E E “bf E < )»?. 

2=0 j — 1 s =0 


(29) 


By subtracting (29) from (28) we obtain Q(r) and Q(r) = |Q(r)|. If Q(r) is 
greater than TO LI, we perform the refinement and use interpolation to estimate 
the values at the new collocation points. 

We study the evolution of the mean and the variance of the solution to 
(22). There is no discontinuity involved in this problem so we can start the 
evolution with only one element. However as time increases the lower ordered 
gPC solution (via Galerkin or collocation approach) starts to deviate from the 
exact solution. One way to keep the error under control is to increase the order 
of the gPC expansion in Galerkin approach or increase the number of collocation 
points in the collocation approach. For the adaptive mesh refinement algorithm 
we have chosen a low order gPC expansion for each element (fewer number of 
collocation points respectively). For the particular example, the result from the 
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Galerkin method is the same as that from collocation method. Here we present 
only results for the collocation approach. Figure 1 shows the curves of the mean 
(left) and variance (right) via various methods. The results from the refinement 
algorithm almost reproduce the exact solution while the global gPC solution 
starts departing from the exact solution as time evolves. We choose the mean 
and variance of the exact solution as references to study the relative error of 
each algorithm. We define 


Error of mean = 


H(u(t;u})) - n(u(t; u>)) 


Error of variance = 


a 2 (u(t ; w)) — cr 2 (u(t ; ui)) 
a 2 (u(t; w)) 


In Figure 2 we present the evolutions of the error for the gPC solution with order 
5, for the adaptive Galerkin approach of order 5 and the adaptive collocation 
approach of order 5 for different values of the accuracy control parameters. 
The maximum relative errors for the mean and the variance resulting from the 
adaptive mesh refinement with different orders are presented in Table 1. For 
reasons of comparison we include the relative error from the gPC solution of 
order 5. As expected the higher ordered models require fewer elements while at 
the same time they maintain higher accuracy. The adaptive meshes at t = 10 
with different accuracy control values TOL\ = 10“ 1 and TOL\ = 10 -2 are 
demonstrated in Figure 3. The elements on the left end are smaller than those 
on the right end. This is consistent with the fact that the rate of change of 
u(t; u>) is larger on the left end of the random space (this is because for negative 
values of the parameter the solution blows up exponentially and thus is more 
sensitive to uncertainty). The evolution of the error, of the number of elements 
and of the error for the variance are shown in Figure 4. 




Figure 1: Evolution of mean of «(f; oj) (left.) and evolution of variance of u(t; uj) (right) for the 
simple ODE. 
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Figure 2: Evolution of relative error of mean of u(t; uj) (left) and evolution of relative error of 
variance of u(t] uj) (right) for the simple ODE. 



N 

Error of p(u) 

Error of var(u) 

gPC, p = 5 

1 

3.8e — 3 

l.le — 1 

TOLi = 1CT 1 

p = 5 

15 

7.3e-5 

5.7e — 4 

P= 7 

9 

1.5e — 6 

3.3e — 5 

TOL\ = 1(T 2 

p = 5 

19 

1.0e-5 

8.0e — 5 

P = 7 

11 

3.0e- 7 

5.6e — 6 


Table 1: Maximum of relative error of mean and variance of solution to the simple ODE when 
t e (0,10] 



Figure 3: Result of adaptive mesh refinement for p = 7, TOL\ = 10 1 : length of the elements 
at t = 10 (left) and distribution of the collocation points at t = 10 (right). 
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Figure 4: Evolutions of mean error and variance error compared with adaptive meshes for the 
simple ODE when p = 7, TOL\ = 10 -2 . 


4-2. The Kraichnan-Orszag three-mode system 

Consider the following system obtained by a linear transformation performed 
on the original Kraichnan-Orszag (K-O) three-mode system [16, 18, 21], 

dyi 
dt 
dyi 
dt 
dy 3 
dt 

subject to initial conditions, 

2/1 (0) = 3/i (0;w), y 2 (0) = y 2 (0; w), 2/3(0) =2/ 3 (0; w). (31) 

The solution of this system changes qualitatively depending on the value of the 
initial conditions 2/1 (0) and 2 / 2 ( 0 ). For certain initial conditions the solution 
of the problem is periodic, and the period goes to infinity as 2/1 (0) and 2 / 2 ( 0 ) 
approach to the planes 2/1 = 0 and 1/2 = 0. Thus a discontinuity in the initial 
condition dependence exists when the initial conditions are allowed to be on both 
sides of these two planes. It was shown that a global Wiener-Hermite expansion 
does not faithfully represent the dynamics of the system when the random inputs 
are Gaussian variables in [26]. Mesh refinement algorithms [16, 18] and the 


= 2/12/3, 

= -2/22/3, ( 30 ) 

2 1 2 

= -y 1 + 2 / 3 , 
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adaptive hierarchical sparse grid collocation [21] can efficiently quantify the 
uncertainty of the system when the random inputs are uniform random variables 
and are in a range which includes the discontinuity value. We consider the case 
with one random input, two random inputs and three random inputs respectively 
and the random initial conditions coming from a uniform distribution. We 
present results for both the Galerkin and collocation approaches. The associated 
system of ordinary differential equations is solved using the standard fourth- 
order Runge-Kutta scheme with a timestep At = 0.01. 

4-2.1. One-dimensional random input 
We choose the initial conditions 

S/i(0;w) = l, j/2(0;w) = 0.1£(a;), y 3 (0;cd)=0, (32) 

where £ ~ U[— 1,1]. In this case the discontinuity point z /2 = 0 is in the random 
input space. 




Figure 5: Evolution of the variance of y\ (left) and evolution of the variance of y 2 (right) for 
the Kraichnan-Orszag three-mode system with ID random initial condition. 

We study the evolution of each variable jji. i = 1,2,3 on the time interval 
[0, 30] (short time behavior). Figure 5 presents the variance evolution of y\ and 
that of y 2 estimated by the quasi-Monte Carlo Sobol (MC-SOBOL) sequence 
with 1,000,000 samples, AMR-GA with polynomial bases up to order 5 under 
various accuracy control values TOL\ and global gPC Galerkin method with 
order 30. Failure to capture the properties via the global gPC expansion after 
a short time (t >= 8) is also shown in Figure 5. Meanwhile, the results from 
AMR-GA converge. If we keep the order of the expansion fixed and make the 
tolerance TOLi stricter (smaller), more elements are desired via AMR-GA and 
higher accuracy is achieved. 

We consider the solution given by MC-SOBOL with 1,000,000 samples as 
the reference (ref) to study the error of AMR-CO, the global collocation method 
(CO), brute force Monte Carlo (MC) and MC-SOBOL. Results from AMR-CO 
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AMR-CO 

CO 

MC 

MC-SOBOL 



No. of points 

Error 

Error 

Error 

Error 

TOL i 

= 1CT 3 

160 

4.6e — 2 

1.7e — 1 

3.0e — 1 

2.7e — 1 

TOL 1 

= 10 -4 

260 

4.1e — 3 

1.3e — 1 

2.4e - 1 

1.4e — 1 

TOLi 

= icr 5 

320 

2.8e — 4 

l.Oe — 1 

9e — 2 

l.le — 1 


Table 2: Maximum relative errors of the variance when t E (0, 30] for the Kraichnan-Orszag 
three-mode system with ID random input via different method 



N 

Error 

N 

Error 

N 

Error 

TOL\ 

10~ 3 

10 -5 

io~ 7 

P = 7 

30 

1.7e — 1 

44 

1.8e — 4 

86 

5.0e — 6 

p = 9 

20 

9.7e — 2 

34 

2.1e — 4 

62 

6.8e — 7 

p = 11 

30 

l.le- 1 

34 

3.7e — 4 

52 

1.7e — 6 


Table 3: Maximum relative errors of the variance when t E (0, 30] for the Kraichnan-Orszag 
three-mode system with ID random input via AMR-GA 



N 

Error 

N 

Error 

N 

Error 

TO Lx 

10” 3 

10“ 5 

10~ 7 

P= 7 

182 

3.8e — 2 

352 

2.2e — 4 

688 

5.1e — 6 

p = 9 

160 

4.6e — 2 

320 

2.8e — 4 

640 

9.9e — 7 

p =11 

216 

8.4e - 2 

312 

2.4e - 4 

624 

2.0e — 6 


Table 4: Maximum relative errors of the variance when f G (0, 30] for the Kraichnan-Orszag 
three-mode system with ID random input via AMR-CO 
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(a) (b) 


Figure 6: (a) Adaptive meshes for the Kraichnan-Orszag three-mode system with ID random 
input when p = 7 and TOL\ = 10“ 3 via AMR-GA; (b) zoom-in mesh of (a) near £ = 0. 





TOL^IoAn = 192 | 




(a) (b) 


Figure 7: (a) Distribution of collocation points when solving the Kraichnan-Orszag three¬ 
mode system with ID random input via AMR-CO with p = 7 and TOL i = 10 3 ; (b) zoom-in 
mesh of (a) near £ = 0. 


with p = 9 and different accuracy control values(TOLi) are presented in Table 
2. The error is defined as: 


Error = 


sup max 
te(0,30] i=1 > 2 > 3 


|var(y,)(t) - var ref (^)(£)| 
var re f(yi)(t) 


(33) 


In the first column of Table 2 we present the number of collocation points and 
the maximum of the relative errors of the variance of y±, j /2 and ys when the 
system evolves to t, = 30 via AMR-CO with full model orders 9 under different 
TOLi. In the implementation of the global collocation method we use Legendre 
collocation points with the same number of points as in AMR-CO. For MC and 
MC-SOBOL, the sample size equals the number of the collocation points in 
AMR-CO. Errors are shown in column two, three and four respectively. We 
note a much higher convergence rate for AMR-CO than for the other methods. 

Secondly, since for this example the mesh refinement method appears to 
have a higher convergence rate we choose the result from AMR-CO with order 
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19 and very strict accuracy control TOL\ = 10~ 9 as the reference to study 
the relative error of both adaptive mesh refinement methods (AMR-GA and 
AMR-CO) with set of different orders and TOL\. Table 3 shows the results of 
AMR-GA while Table 4 presents the results of AMR-CO. For a fixed value of 
the tolerance TOL \, higher order models require fewer elements in AMR-GA. 
However if we check carefully the number of variables we resolve in AMR-GA 
at t = 30 for each order, the magnitude of numbers are comparable for fixed 
TOL i. For example, when TOL\ = 10~ 7 , we solve 688 variables at t = 30 for 
order 7 model, while that number becomes 620 for order 9 model and 624 for 
order 11 model respectively, and the relative errors of all three orders are at 
level 10~ 6 . Furthermore, the same trend appears in the results of AMR-CO. 

Details of the adaptive meshes resulting from our AMR-GA and AMR-CO 
algorithms around £ = 0 are presented in Figures 6 and 7. The finest meshes 
are around the discontinuity of the random space. This shows that our mesh 
refinement criterion locates accurately the discontinuity. Furthermore, when 
TOL\ is extremely small, the meshes exhibit the pattern that the closer the 
element is to £ = 0, the smaller the element in AMR-GA and the more the 
collocation points in AMR-CO. Meanwhile the meshes and points are symmetric 
with respect to £ = 0 as they should be according to the symmetry of the K-0 
system. 

4-2.2. Two-dimensional random input 

We study the K-0 system with initial conditions involving two independent 
random inputs 

2/i(0;w) = l, y 2 (0;w) = 0.1£i(w), y 3 (0; w) = &M, (34) 

where and £2 are independent uniform random variables on [—1,1]. In Figure 
8, we plot the evolution of the variance of each random output yi,i = 1,2,3 
subject to a 2D random input obtained by AMR-GA with order 7 for the full 
model with different values of TOLi. We also show the mesh of the random 
space at time t = 10 for 210 elements generated by AMR-GA with order 7 for 
the full model and, TOL\ = 10 -6 and TOL 2 = 0.1. The smallest elements are 
around the discontinuity y 2 = 0 and the results are more sensitive to £1 because 
of the discontinuity introduced by £ 1 . 

To study the convergence of AMR-CO compared with CO, MC and MC- 
SOBOL, we take the result from MC-SOBOL with 1,000,000 samples as the 
reference. To implement the direct collocation method, we choose the tensor grid 
with \N 2 ] if [IVis even, otherwise we choose [IV 2 ] -f 1 Legendre quadrature 
points in each dimension. As for MC and MC-SOBOL, we take the sample size 
to be the number of collocation points in AMR-CO. Table 5 demonstrates the 
relative error with respect to the reference. The results show that AMR-CO 
converges much faster than the other three methods. 

The result from AMR-CO with order 10 for the full model, TOL\ = 10~ 9 
and TOL 2 = 0.1 is selected to be the reference to derive the relative errors of 
low ordered models. Table 6 presents the relative errors of variance by AMR- 
GA of different sets of full model orders and different levels of accuracy control, 
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while Table 7 presents the relative errors of variance by AMR-CO. We observe 
similar trends as in the ID case, namely that more accurate models require 
fewer elements in AMR-GA if the accuracy tolerance is kept fixed. Also, that 
stricter accuracy control requires more elements if the order of the models is 
fixed. To achieve the same level of relative errors, the number of the elements 
increases faster in the 2D case than the ID case. However, in AMR-CO with 
the same accuracy control, the number of required collocation points are in the 
same order of magnitude for different full model orders. Also, the relative errors 
are in the same order of magnitude. 


AMR-CO 

CO 

MC 

MC-SOBOL 



No. of points 

Error 

Error 

Error 

Error 

TOLi = 

10~ 2 

576 

1.8e — 1 

3.7e — 1 

1.5e — 1 

2.5e — 1 

TOLi = 

10 -3 

1944 

3.7e — 2 

1.7e — 1 

l.Oe — 1 

2.1e — 1 

TOLi = 

10 -4 

4896 

3.1e — 3 

8.5e — 2 

6.2e — 2 

9.3e — 2 


Table 5: Maximum relative errors of the variance when t £ (0,10] for the Kraichnan-Orszag 
three-mode system with 2D random input via different method 



N 

Error 

N 

Error 

N 

Error 

TOL i 

10" 3 

10 -5 

10~ 7 

p = 5 

34 

2.8e — 2 

222 

2.6e — 3 

424 

1.7e — 4 

P = 7 

32 

2.7e — 2 

152 

7.7e - 4 

310 

4.7e — 6 


Table 6: Maximum relative errors of the variance when t £ (0,10] for the Kraichnan-Orszag 
three-mode system with 2D random input via Galerkin approach 



N 

Error 

N 

Error 

N 

Error 

TOTi 

10~ 3 

10 -5 

10~ 7 

p = 5 

1944 

3.7e — 2 

10224 

3.2e — 4 

26496 

3.8e — 6 

P = 7 

2048 

7.4e - 2 

9728 

8.9e — 4 

19840 

6.4e — 6 


Table 7: Maximum relative errors of the variance when t £ (0,10] for the Kraichnan-Orszag 
three-mode system with 2D random input via collocation approach 


4-. 2.3. Three-dimensional random input 
The initial conditions in this case are 

J/i (0;w) = £i(w), 2/ 2 (0 ;w) = 6M, 2/3(0; w) = £ 3 (w), (35) 
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( a ) 




(b) 



(c) 


(d) 


Figure 8: The Kraichnan-Orszag three-mode system with 2D random inputs, p = 7, TOL 2 = 
0.1. (a) Evolution of variance of y\\ (b) Evolution of variance of 1 / 2 '- (c) Evolution of variance 
of y ? } : (d)Adaptive meshes for 2D random input when p = 7, TOLi = 10 6 , N = 210. 


where £ 1 , £2 and £3 are independent uniform random variables on [—1,1], In 
this case, discontinuities occur at y\ = 0 and t /2 = 0. Due to the presence 
of discontinuities and moderately high dimensionality, this problem is more 
difficult than the previous ones. Figure 9 shows the evolution of the variance of 
y\ and y 3 obtained by an order 5 full model with different TOL\. The variances 
of y\ and t /2 are the same due to the symmetry of y\ and y 2 in (30) and the 
corresponding initial random input. The results for a global gPC expansion of 
order 9 diverges from the MC-SOBOL results at t ps 3. On the other had, our 
AMR-GA algorithms achives much better accuracy. Similarly as in the previous 
two cases, a comparison of the relative errors of the output variances between 
AMR-CO, global collocation, MC and MC-SOBOL is presented in Table 8 . (the 
MC-SOBOL result with 1,000,000 samples is taken as the reference). Again, the 
mesh refinement algorithm shows a better convergence than the other methods, 
however the advantage is not as big as before. 

We choose results from AMR-CO with p = 10, TOL 1 = 10 ~ 7 as the reference 
and show the results of AMR-GA with different set of orders and TOL\ in Table 
9 and that of AMR-CO in Table 10. As we can see, the number of the elements in 
AMR-GA as well as the collocation points in AMR-CO grows dramatically fast 
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in order to gain sufficient accuracy compared to the ID and 2D cases. Therefore, 
both of our adaptive mesh refinement algorithms are still dimension-dependent 
algorithms and many more elements or collocations points are required to resolve 
the discontinuity. 




Figure 9: Evolution of the variance of yi(left) and evolution of the variance of ys (right) for 
the Kraichnan-Orszag three-mode system with 3D random inputs. 


AMR-CO 

CO 

MC 

MC-SOBOL 



No. of points 

Error 

Error 

Error 

Error 

TOL 1 

= 1CT 2 

2784 

3.4e - 2 

8 . 1 e — 2 

4.7e — 2 

6.4e - 2 

TOL 1 

= icr 3 

4176 

2.4e - 2 

7.0e - 2 

3.7e — 2 

8.5e — 2 

TOLi 

= icr 4 

23664 

3.4e - 3 

3.8e - 2 

1 . 6 e — 2 

2 . 8 e — 2 


Table 8: Maximum relative errors of the variance when t € (0, 6] for the Kraichnan-Orszag 
three-mode system with 3D random input via different method 


4-3. Kuramoto-Sivashinsky equation 

In this section, we study the normalized Kuramoto-Sivashinsky (K-S) equa¬ 
tion with periodic boundary conditions in the interval [0, 27r]. Let the damping 
parameter to the original value derived by Sivashinsky, v = 4, and introduce 

• • — n r 2 ^ 

the bifurcation parameter a = 4 L = j -. The equation can be restated as 


u t = -4 u x 


: + jW 2 


0 < x < 27r, 


u(x + 27t, t) = u(x, t), u(x, 0) = uq(x). 


(36) 


The K-S equation solution exhibits a broad range of behaviors as a varies. In 
particular, as a. is varied the solution can exhibit bifurcations, thus a small 
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N 

Error 

N 

Error 

N 

Error 

TOL\ 

10 -2 

10 -3 

10" 4 

p = 4 

32 

8 .6e — 2 

48 

3.0e - 2 

312 

2.7e — 3 

p = 5 

32 

7.1e — 2 

40 

5.3e — 2 

280 

3.0e — 3 

p = 6 

16 

9.9e — 2 

32 

2.9e — 2 

112 

1.7e — 3 


Table 9: Maximum relative errors of the variance when t £ (0, 6] for the Kraichnan-Orszag 
three-mode system with 3D random input via Galerkin approach 



N 

Error 

N 

Error 

N 

Error 

TOLi 

10” 3 

10 -5 

io- 7 

p = 4 

2784 

3.4e - 2 

4176 

2.4e — 2 

23664 

3.4e — 3 

p = 5 

2160 

5.9e — 2 

5400 

7.9e - 3 

29160 

3.1e — 3 

p = 6 

3312 

1.3e — 1 

6624 

2 .1e — 2 

21528 

3.2e — 3 


Table 10: Maximum relative errors of the variance when t £ (0, 6] for the Kraichnan-Orszag 
three-mode system with 3D random input via collocation approach 


change in a may result in a dramatic change of the solution. We examine the 
mean of the solution when the parameter a is in an interval which includes 
several bifurcation values. High precision is required for the numerical simula¬ 
tion due to the extreme sensitivity of the solution to the stepsize used [27, 28]. 
We use high precision pseudospectral approximations for the spatial derivatives, 
and semi implicit time differencing with fixed time step for the temporal deriva¬ 
tive. We only present results for the adaptive mesh refinement algorithm with 
collocation (AMR-CO). 

We restrict attention to the interval a £ [13,17] in which the solution can be 
unimodal or bimodal depending on the value of a. We obtain the solution for 
the initial condition uq(x) = 2.6680cos(a:) + 0.1979cos(2:r)+ 0.0094cos(3a;). For 
a in the unimodal regime, the solution quickly converges to the unimodal steady 
solution. Figure 10(a) shows the unimodal solution when a = 13 while Figure 
10(b) presents the periodical bimodal solution when a = 17. The solution to 
K-S equation is very sensitive with respect to the parameter a, and we need 
to start with an adequate number of elements in random space. We start with 
32 initial elements in each refinement implementation. Since the expression 
of the RHS in (36) is complicated, only the decoupled AMR-CO is simulated. 
Result from MC-SOBOL simulation with 1,000,000 samples is considered as 
the reference. 

Figure 11 shows the evolution of the mean of K-S solution vis AMR-CO with 
d = ll,TOLi = 0.1 compared with that of the reference. Figure 12 shows the 
comparison of the evolution of the variance of K-S solution between AMR-CO 
and the reference MC-SOBOL. We can see that although the behavior of the 
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solution is very complicated, the AMR-CO still captures the major properties 
of the solution. The evolution variances at x = 5.1296 by different orders 
compared with the reference is shown in Figure 13(a). In this case, a more 
actuate model(e? = 11) results in a better simulation. The distributions and 
length of the elements of the parameter space are presented in Figure 13(b). 
Since when d = 11, there are 12 collocation points in each clement, in this case 
more smaller elements as well as more points are distributed in the bimodal 
solution regime. The final mesh has a lot of elements which demonstrate extreme 
sensitivity of the solution to the bifurcation parameter. 



(a) 


(b) 


Figure 10: Solutions to K-S equation with bifurcation parameter a = 13(left) and a = 
17(right) when t E [0,10] 



(a) 


(b) 


Figure 11: Mean of solution to K-S equation with bifurcation parameter a ~ C/ [13, 17] when 
t G [0,10] via AMR-CO with d = 11, TOL\ = O.l(left) and MC-SOBOL with 1,000,000 
samples (right) 


4-4- Refinement in physical space 

Even though we have presented the formulation of our mesh refinement algo¬ 
rithm for random space, the same framework can accommodate also refinement 
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(a) (b) 

Figure 12: Variance of solution to K-S equation with bifurcation parameter a ~ U[ 13,17] 
when t G [0,10] via AMR-CO with d = 11, TOL\ — 0.1 (left) and MC-SOBOL with 1, 000, 000 
samples (right) 




(a) (b) 

Figure 13: Variance of solution to K-S equation at x = 5.1296 (left), Elements of parameter 
space at t = 10 via AMR-CO with d = 11, TOL\ = O.l(right) 


in physical space. We will present here results for the location of shocks for the 
ID Burgers equation (more elaborate examples will be presented in a forthcom¬ 
ing publication). The ID inviscid Burgers equation reads 

u t + uu x = 0, x e (—1,1), (37) 

with initial condition u(x, 0) = uq(x) and boundary condition u(—1, t) = it( 1, t) = 
0. We use the initial condition: 

u(x, 0) = sin(27ra;), x € (— 1,1). (38) 

The solution develops two shocks in finite time. We employ a collocation method 
and solve the resulting system of ODEs by a fourth-order Runge-Kutta scheme 
with At = 10~ 5 .. Figure 14(a) shows the solution of (37) at t = 0.1592(ss l/-7r) 
using a global collocation method with 256 Gauss-Legendre collocation points. 
Because the distribution of the collocation points are concentrated near the 
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boundaries of the domain as seen in top figure in Figure 15, the solution shows 
oscillations in the areas where the gradient of the solution is large. Figure 14(b) 
shows the solution at t = 0.1592 using the adaptive mesh refinement collocation 
algorithm. To resolve the solution we start with 8 elements and 6 collocation 
points in each element. As shown in Figure 15 the adaptive mesh refinement 
algorithm distributes more elements where the slope is large, which means more 
points near the location of the shocks. As a result the oscillations are eliminated. 
Note that the assignment of more points where a better resolution is needed 
stabilizes the numerical scheme. Even with large time step At = 10~ 2 , the mesh 
refinement algorithm is still stable, while direct collocation with 256 Gauss- 
Legendre point becomes unstable and oscillates a lot. 



(a) 


.TOL,.10-*.N-48 



(b) 


Figure 14: Solution of burgers equation by 256 Legendre collocation points(left) and mesh 
refinement algorithm with p = 5, TOL\ = 10 -2 ,10 -4 (right), here At = 10 -5 


_i_i_ 


_i_i— 


i-1 


Figure 15: Distributions of 256 Legendre collocation points (top) and collocation point from 
mesh refinement algorithm with p = 5, TOL\ = 10 2 ,10“ 4 (bottom) 


5. Discussion and future work 

We have presented a novel method for adaptive mesh refinement in the con¬ 
text of uncertainty quantification when a spectral decomposition of the solution 


26 






















is possible. The layering structure of the spectral representation of the solution 
can capture the transfer of activity across scales and thus can be utilized to 
detect when and where higher resolution is needed. This idea is inspired by [24] 
where the reduced model of the full system can be constructed. In this case 
the memory term of the reduced model captures the energy transfer between 
scales and provides a reliable indicator for the mesh refinement process. The 
proposed approach was implemented in the context of multi-element generalized 
polynomial chaos expansions for both Galerkin and collocation methods. 

The same mesh refinement framework can be immediately extended to cases 
requiring refinement in physical space. However, unlike the random space where 
each element can be considered independently, when differentiating or integrat¬ 
ing in physical space all elements are coupled due to the need of regularity of 
the solution across the elements. In the current work we have only included 
results from the application of the mesh refinement algorithm to the accurate 
location of shocks for the ID Burgers equation. More elaborate examples (for 
more spatial dimensions) will be presented in a forthcoming publication. 
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